
#include "wrapper.h"

#include "libmesh.h"
#include "utility.h"
#include "auto_ptr.h"
#include "petsc_vector.h"
#include "petsc_matrix.h"
#include "parameters.h"

//#include "petsc_vec_utils.h"
#include "neutronics_system.h"

namespace Wrapper
{
	bool jac_context = false ;
	bool prec_context = false ;
	bool prec_prec_context = false ;
	bool prec_created = false ;
	double prec_atol=1e-8, prec_rtol=1e-1 ;

	bool initialized = false ;
	AutoPtr<PetscVector<Real> > work1, work2, work3 ;

	orthonorm_function fn_ref = NULL ;

#undef __FUNCT__
#define __FUNCT__ "initialize_wrapper"
	void initialize_wrapper (const PetscVector<Real>& x)
	{
		if (!initialized)
		{
			PetscErrorCode ierr=0;

			work1 = AutoPtr<PetscVector<Real> >(dynamic_cast<PetscVector<Real>*>(x.clone().release())) ;
			work2 = AutoPtr<PetscVector<Real> >(dynamic_cast<PetscVector<Real>*>(x.clone().release())) ;
			work3 = AutoPtr<PetscVector<Real> >(dynamic_cast<PetscVector<Real>*>(x.clone().release())) ;
			//work1 = x.clone() ;
			//work2 = x.clone() ;
			//work3 = x.clone() ;

			initialized = true ;
		}
	}

#undef __FUNCT__
#define __FUNCT__ "clear_wrapper"
	void clear_wrapper ()
	{
		if (initialized)
		{
			work1.reset() ;
			work2.reset() ;
			work3.reset() ;

			initialized = false ;
		}
	}

	// solve Loss operator by repeatedly applying the block-gauss-seidel preconditioner
#undef __FUNCT__
#define __FUNCT__ "__solve_loss_matrix_exactly__"
	PetscErrorCode __solve_loss_matrix_exactly__ (NeutronicsSystem *system, Vec b, Vec x, int max_it=1e3, Real tolerance=1e-12, bool errorcheck=true)  
	{
		PetscErrorCode ierr=0;

		PetscMatrix<Real>* Loss = system->get_loss_matrix() ;

		Vec res, work ;
		ierr = VecDuplicate(b, &res) ;CHKERRQ(ierr);
		ierr = VecDuplicate(b, &work) ;CHKERRQ(ierr);

		// ierr = VecCopy(b, x) ;CHKERRQ(ierr);

		ierr = VecSet(x, 0.0) ;CHKERRQ(ierr);

		Real omega_max = 1.3 ;
		PetscBool monitor_flag = PETSC_FALSE, const_omega = PETSC_FALSE ;

		ierr = PetscOptionsGetReal("lossoptr_", "-ksp_atol",&tolerance,PETSC_NULL);
		ierr = PetscOptionsGetInt("lossoptr_", "-ksp_max_it",&max_it,PETSC_NULL);
		ierr = PetscOptionsGetReal("lossoptr_", "-ksp_omega",&omega_max,PETSC_NULL);
		ierr = PetscOptionsHasName("lossoptr_", "-const_omega",&const_omega);
		ierr = PetscOptionsHasName("lossoptr_", "-monitor",&monitor_flag);

		Real errnrm0 = 0.0, error = 1e2 ;
		Real omega = 1.0 ;

		// find initial residual
		ierr = MatMult(Loss->mat(), x, res) ;CHKERRQ(ierr);
		// res = b - L*x
		ierr = VecAYPX(res, -1.0, b) ;CHKERRQ(ierr);

		ierr = VecNorm(res, NORM_2, &errnrm0) ;CHKERRQ(ierr);
		error = errnrm0 ;

		if (monitor_flag)
		{
			std::cout << "Loss operator solve: tolerance = " << tolerance << ", ||residual_0|| = " << errnrm0 ;
		}

		Real lambda_estimate ;
		Real spectral_radius = 1.0, preverror = errnrm0, omega_estimate = 1.0, nrmx = 0 ;
		int iter = 0;
		while ( error > tolerance*errnrm0 && iter < max_it)
		{
			ierr = precondition_loss_matrix(system, res, work) ;CHKERRQ(ierr);

			ierr = VecAXPY(x, omega, work) ;CHKERRQ(ierr);

			// calculate the new residual
			ierr = MatMult(Loss->mat(), x, res) ;CHKERRQ(ierr);

			// lambda_estimate = || L * x ||
			ierr = VecNorm(res, NORM_2, &lambda_estimate) ;CHKERRQ(ierr);
			ierr = VecNorm(x, NORM_2, &nrmx) ;CHKERRQ(ierr);
			// lambda_estimate = || L * x || / || x || = crude estimate of maximum eigen value
			lambda_estimate = lambda_estimate/nrmx ;

			// res = b - L*x
			ierr = VecAYPX(res, -1.0, b) ;CHKERRQ(ierr);

			ierr = VecNorm(res, NORM_2, &error) ;CHKERRQ(ierr);

			spectral_radius = error/preverror ;

			if (spectral_radius > 1)
			{
				omega_estimate = 0.9 ;
			}
			else
			{
				// std::cout << "|convergence rate| = " << spectral_radius << ", |omega_est| = " << omega_estimate 
				// << ", |omega_val| = " << 2.0/(1.0+sqrt(fabs(1.0-spectral_radius))) ;

				// http://www.jstor.org/stable/2153527?seq=3
				// And Bahi J., et al. Parallel iterative algorithms.. from sequential to grid computing 
				// (ISBN 9781584888086)(Taylor, 2008) : Page 38
				// spectral_radius : 1.0 / ( lowest eigenvalue estimate )
				double omega_estimate1 = 2.0/(1.0+sqrt(fabs(1.0-spectral_radius*spectral_radius))) ;

				// double lambda_dominance = log(spectral_radius) ;
				double omega_estimate2 = 2.0/(1.0+sqrt(fabs(1.0-spectral_radius*lambda_estimate))) ;

				omega_estimate = (omega_estimate1 > omega_estimate2 ? omega_estimate1 : omega_estimate2) ;

				// std::cout << "Iteration [" << iter << "]\t||residual|| = " << error 
				// << ",  |spectral_radius_est| = " << 1.0/spectral_radius 
				// << ", |omega_est_1| = " << omega_estimate1 << ", |omega_est_2| = " << omega_estimate2 
				// << ", |omega_est| = " << omega_estimate << ", |omega_max| = " << omega_max ;
			}

			if (const_omega)
			{
				omega_estimate = omega_max ;
			}

			iter++ ;
			preverror = error ;

			if (monitor_flag)
			{
				std::cout << "Iteration [" << iter << "]\t||residual|| = " << error 
					<< ",  |spectral_radius_est| = " << 1.0/spectral_radius 
					<< ", |omega_est| = " << omega_estimate << ", |omega_max| = " << omega_max ;
			}

			if (omega_estimate > omega_max)
				omega_estimate = omega_max ;

			omega = omega_estimate ;
		}

		if (error < tolerance*errnrm0 && iter < max_it && monitor_flag)
		{
			std::cout << "--Iteration converged--" ;
		}

		if (errorcheck && monitor_flag)
		{
			if (iter > max_it)
			{
				PetscPrintf(PETSC_COMM_WORLD, "--Loss matrix inversion through SOR Iteration: Exceeded maximum iterations.--") ; 
				std::cin.get() ;
			}
			else if(error > tolerance*errnrm0)
			{
				PetscPrintf(PETSC_COMM_WORLD, "--Loss matrix inversion: Not converged: ||residual_final|| = %G after %D iterations", error, iter) ;
				std::cin.get() ;
			}
		}

		ierr = VecDestroy(&res) ;CHKERRQ(ierr);
		ierr = VecDestroy(&work) ;CHKERRQ(ierr);

		return ierr ;
	}

	// specifically a method to copy [phi; k] vector to a [phi] only vector
	// and vice versa : Very useful for simple eigen-problem calculations
#undef __FUNCT__
#define __FUNCT__ "petsc_copy_vector"
	PetscErrorCode petsc_copy_vector(Vec xin, Vec xout, bool is_input_smaller)
	{
		PetscInt ldim, low_in, high_in ;
		PetscInt low_out, high_out ;
		VecGetOwnershipRange(xin,&low_in,&high_in);
		VecGetOwnershipRange(xout,&low_out,&high_out);

		libmesh_assert(low_in == low_out) ;

		PetscErrorCode ierr;

		if (is_input_smaller)
			ierr = VecGetLocalSize(xin,&ldim);
		else
			ierr = VecGetLocalSize(xout,&ldim);
		CHKERRQ(ierr);

		PetscScalar *xarr, *yarr ;

		ierr = VecGetArray(xin, &xarr) ;CHKERRQ(ierr);
		ierr = VecGetArray(xout, &yarr) ;CHKERRQ(ierr);

		for (int i=0; i<ldim; i++)
		{
			yarr[i] = xarr[i] ;
		}

		ierr = VecRestoreArray(xin, &xarr) ;CHKERRQ(ierr);
		ierr = VecRestoreArray(xout, &yarr) ;CHKERRQ(ierr);

		return ierr ;
	}


#undef __FUNCT__
#define __FUNCT__ "simple_eigen_operator_private"
	PetscErrorCode simple_eigen_operator_private(NeutronicsSystem *system, Vec x, Vec y) 
	{
		PetscErrorCode ierr = 0;

		// level 3 call -> use work3

		PetscMatrix<Real>* Gain = system->get_gain_matrix() ;

		// pT1 = P * x ; 
		ierr = MatMult(Gain->mat(), x, work3->vec()) ;CHKERRQ(ierr);

		if (!prec_context && !jac_context && !prec_prec_context)
		{
			PUSH_LOG("simple_eigen_operator", "Eigen Context");

			PetscBool flg ;
			ierr = PetscOptionsHasName(PETSC_NULL,"-exact_loss",&flg); ;

			if (flg)
			{
				KSP& loss_solver = system->loss_matrix_solver() ;
				ierr = KSPSolve(loss_solver, work3->vec(), y) ;CHKERRQ(ierr);
			}
			else
			{
				ierr = __solve_loss_matrix_exactly__ (system, work3->vec(), y)  ;CHKERRQ(ierr);
			}

			// VecView(y, 	PETSC_VIEWER_STDOUT_WORLD ) ;
			// std::cin.get() ;

			POP_LOG("simple_eigen_operator", "Eigen Context");
		}
		else
		{
			// it is either jac_context or prec_context
			if (prec_prec_context)
			{
				PUSH_LOG("preconditioners_preconditioner_operator", "Eigen Context");

				unsigned int prec_method = 3 ;

				if (prec_method < 3)
				{
					// PetscMatrix<Real>* Gain = system->get_gain_matrix() ;
					PetscMatrix<Real>* Loss = system->get_loss_matrix() ;

					if (prec_method == 1)
					{
						// pT1 = rowsum(L) ; lumped loss matrix
						Vec work_1 ;
						ierr = VecDuplicate(work3->vec(), &work_1) ;CHKERRQ(ierr);

						VecSet(work_1, 1.0) ; 
						ierr = MatMult(Loss->mat(), work_1, y) ;CHKERRQ(ierr);
						ierr = VecDestroy(&work_1) ;CHKERRQ(ierr);
					}
					else
					{
						// pT1 = rowsum(L) ; diagonal loss matrix = point jacobi
						MatGetDiagonal(Loss->mat(), y) ; 
					}

					ierr = VecPointwiseDivide(y, work3->vec(), y) ;CHKERRQ(ierr);
				}
				else
				{
					ierr = precondition_loss_matrix(system, work3->vec(), y) ;CHKERRQ(ierr);
				}

				POP_LOG("preconditioners_preconditioner_operator", "Eigen Context");
			}
			else if (prec_context)
			{
				PUSH_LOG("simple_preconditioner_operator", "Eigen Context");

				int max_prec_its = 1 ;
				ierr = PetscOptionsGetInt("lossoptr_", "-prec_max_it",&max_prec_its,PETSC_NULL);

				ierr = __solve_loss_matrix_exactly__ (system, work3->vec(), y, max_prec_its, 1e-5, false)  ;CHKERRQ(ierr);

				// ierr = precondition_loss_matrix(system, work3->vec(), y) ;CHKERRQ(ierr);

				POP_LOG("simple_preconditioner_operator", "Eigen Context");
			}
			else
			{
				PUSH_LOG("simple_jacobian_operator", "Eigen Context");

				int max_jac_its = 4 ;
				ierr = PetscOptionsGetInt("lossoptr_", "-jac_max_it",&max_jac_its,PETSC_NULL);

				ierr = __solve_loss_matrix_exactly__ (system, work3->vec(), y, max_jac_its, 1e-10, false)  ;CHKERRQ(ierr);

				// ierr = precondition_loss_matrix(system, work3->vec(), y) ;CHKERRQ(ierr);

				POP_LOG("simple_jacobian_operator", "Eigen Context");
			}
		}

		// apply the shift on the operator if any
		if (fabs(system->eigen_shift()) > 1e-13)
		{
			ierr = VecAXPY(y, system->eigen_shift(), x) ;CHKERRQ(ierr);
		}

		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		PetscFunctionReturn(ierr);
	}

	// y = inv(L)*P*x; Mat 'A' represents the operator = inv(L)*P
#undef __FUNCT__
#define __FUNCT__ "simple_eigen_operator"
	PetscErrorCode simple_eigen_operator(Mat A, Vec x, Vec y)
	{
		void           *ctx;
		NeutronicsSystem *system ;
		PetscErrorCode ierr = 0;

		// A = operator matrix
		ierr = MatShellGetContext(A, &ctx); CHKERRQ(ierr);
		system = (NeutronicsSystem *)ctx ;

		libmesh_assert(system != NULL) ;

		ierr = simple_eigen_operator_private(system, x, y) ;	CHKERRQ(ierr);

		PetscFunctionReturn(ierr);
	}

	// y = inv(L)*P*x; Mat 'A' represents the operator = inv(L)*P
#undef __FUNCT__
#define __FUNCT__ "simple_eigen_preconditioner_solve_operator"
	PetscErrorCode simple_eigen_preconditioner_solve_operator(PC pc, Vec x, Vec y)
	{
		PetscErrorCode ierr = 0;

		void *ctx ;
		ierr = PCShellGetContext(pc, &ctx);CHKERRQ(ierr);

		NeutronicsSystem* system = (NeutronicsSystem *)ctx ;

		libmesh_assert(system != NULL) ;

		prec_context = true ;

		Mat prec_matrix ;

		int numdofs, locdofs;
		ierr = VecGetSize(x, &numdofs); CHKERRQ(ierr);
		ierr = VecGetLocalSize(x, &locdofs); CHKERRQ(ierr);

		ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)system, &prec_matrix);CHKERRQ(ierr);
		ierr = MatSetFromOptions(prec_matrix) ;CHKERRQ(ierr);

		ierr = MatShellSetOperation(prec_matrix, MATOP_MULT, (void(*)())Wrapper::simple_eigen_operator);CHKERRQ(ierr);

		KSP ksp_prec ;

		ierr = KSPCreate(PETSC_COMM_WORLD, &ksp_prec) ;CHKERRQ(ierr);
		ierr = KSPSetType(ksp_prec, KSPGMRES) ;
		KSPAppendOptionsPrefix(ksp_prec,"st_prec_");
		ierr = KSPSetOperators(ksp_prec, prec_matrix, prec_matrix, SAME_NONZERO_PATTERN) ; CHKERRQ(ierr);

		PC pc_prec ;
		ierr = KSPGetPC(ksp_prec, &pc_prec) ;CHKERRQ(ierr);
		ierr = PCSetType(pc_prec, PCNONE) ; CHKERRQ(ierr);

		ierr = KSPSetFromOptions(ksp_prec) ;CHKERRQ(ierr);

		// VecView(x, 	PETSC_VIEWER_STDOUT_WORLD ) ;
		// std::cin.get() ;

		// std::cout << "now printing residual "  ; 

		ierr = KSPSolve(ksp_prec, x, y) ; CHKERRQ(ierr);

		// VecView(y, 	PETSC_VIEWER_STDOUT_WORLD ) ;
		// std::cin.get() ;

		ierr = KSPDestroy(&ksp_prec) ;	CHKERRQ(ierr);
		ierr = MatDestroy(&prec_matrix) ;	CHKERRQ(ierr);

		prec_context = false ;

		PetscFunctionReturn(ierr);
	}

	// y = lambda*x - inv(L)*P*x
#undef __FUNCT__
#define __FUNCT__ "simple_eigen_residual"
	PetscErrorCode simple_eigen_residual(NeutronicsSystem *system, Vec x, Vec y)
	{
		PetscErrorCode ierr = 0;

		int n, nloc ;
		VecGetSize(y, &n) ;

		ierr = petsc_copy_vector(x, work2->vec(), false) ;	CHKERRQ(ierr);

		double vecnrm = work2->l2_norm() ;
		vecnrm *= vecnrm ;

		ierr = simple_eigen_operator_private(system, work2->vec(), work1->vec()) ;CHKERRQ(ierr);

		Real eigval = 1.0 ;
		//		if(KARMA::proc_id() == KARMA::num_procs()-1)
		{
			PetscScalar *xarray ;

			VecGetLocalSize(x, &nloc) ;

			VecGetArray(x, &xarray);

			eigval = xarray[nloc-1] ;

			VecRestoreArray(x,&xarray);
		}

		// res = inv(L)*P*x - lambda*x
		ierr = VecAXPY(work1->vec(), -eigval, work2->vec()) ;	CHKERRQ(ierr);

		ierr = petsc_copy_vector(work1->vec(), y, true) ;	CHKERRQ(ierr);

		// set residual = (0.5-0.5*x*xT) for the eigen value equation
		ierr = VecSetValue(y, n-1, (0.5*vecnrm-0.5), INSERT_VALUES) ;CHKERRQ(ierr);

		// std::cout << "keff = " << eigval << " and neutronics system keff = " << system->Keff() 
		// << " and eigenvalue residual = " << (0.5*vecnrm-0.5) ;

		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		PetscFunctionReturn(ierr);
	}


	// y = L*lambda*x - P*x
#undef __FUNCT__
#define __FUNCT__ "generalized_eigen_residual"
	PetscErrorCode generalized_eigen_residual(NeutronicsSystem *system, Vec x, Vec y, double keig)
	{
		PetscErrorCode ierr = 0;

		unsigned int nloc ;
		VecGetLocalSize(x, (int*)&nloc) ;

		libmesh_assert(nloc == work1->size()) ;

		PetscMatrix<Real>* Gain = system->get_gain_matrix() ;
		PetscMatrix<Real>* Loss = system->get_loss_matrix() ;

		// work1 = P * x ; y = L * x
		ierr = MatMult(Gain->mat(), x, work1->vec()) ;CHKERRQ(ierr);	work1->close() ;

		ierr = MatMult(Loss->mat(), x, y) ;CHKERRQ(ierr);

		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		// res(1:n-1) = P * x - L * lamb * x 
		ierr = VecAYPX(y, -keig, work1->vec()) ;	CHKERRQ(ierr);

		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		PetscFunctionReturn(ierr);
	}

	// y = L*lambda*x - P*x
#undef __FUNCT__
#define __FUNCT__ "generalized_eigen_residual"
	PetscErrorCode generalized_eigen_residual(NeutronicsSystem *system, Vec x, Vec y)
	{
		PetscErrorCode ierr = 0;

		unsigned int nloc ;
		VecGetLocalSize(x, (int*)&nloc) ;

		libmesh_assert(nloc == work1->size()) ;

		PetscMatrix<Real>* Gain = system->get_gain_matrix() ;
		PetscMatrix<Real>* Loss = system->get_loss_matrix() ;

		// work1 = P * x ; pT2 = L * x
		ierr = MatMult(Gain->mat(), x, work1->vec()) ;CHKERRQ(ierr);	work1->close() ;

		ierr = MatMult(Loss->mat(), x, y) ;CHKERRQ(ierr);
		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		// Real keig = system->Keff() ;
		// Use Rayleigh estimate to find Keff
		Real keig = 1.0, production=0.0, loss=0.0 ;
		ierr = VecDot(x, work1->vec(), &production) ; 	CHKERRQ(ierr);
		ierr = VecDot(x, y, &loss) ; 	CHKERRQ(ierr);
		keig = production / loss ;

		// res(1:n-1) = P * x - L * lamb * x 
		ierr = VecAXPY(y, -keig, work1->vec()) ;	CHKERRQ(ierr);

		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		PetscFunctionReturn(ierr);
	}

	// y = L*x; If Mat is shell, return operator application; If Mat is valid, return product
#undef __FUNCT__
#define __FUNCT__ "multiply_loss_matrix"
	PetscErrorCode multiply_loss_matrix(Mat L, Vec x, Vec y)
	{
		void           *ctx;
		NeutronicsSystem *system ;
		PetscErrorCode ierr = 0;

		// L = loss matrix
		ierr = MatShellGetContext(L, &ctx); CHKERRQ(ierr);
		system = (NeutronicsSystem *)ctx ;

		libmesh_assert(system != NULL) ;

		const PetscVector<Real> pX(x) ;
		PetscVector<Real> pY(y) ;

		// direct to the physics solver and make the actual call
		system->loss_matrix_vector(pX, pY) ; 

		PetscFunctionReturn(ierr);
	}

	// y = P*x; If Mat is shell, return operator application; If Mat is valid, return product
#undef __FUNCT__
#define __FUNCT__ "multiply_gain_matrix"
	PetscErrorCode multiply_gain_matrix(Mat P, Vec x, Vec y)
	{
		void           *ctx;
		NeutronicsSystem *system ;
		PetscErrorCode ierr = 0;

		// P = gain matrix
		ierr = MatShellGetContext(P, &ctx); CHKERRQ(ierr);
		system = (NeutronicsSystem *)ctx ;

		libmesh_assert(system != NULL) ;

		const PetscVector<Real> pX(x) ;
		PetscVector<Real> pY(y) ;

		// direct to the physics solver and make the actual call
		system->gain_matrix_vector(pX, pY) ; 

		PetscFunctionReturn(ierr);

	}


#undef __FUNCT__
#define __FUNCT__ "multiply_true_jacobian_matrix"
	PetscErrorCode multiply_true_jacobian_matrix(Mat J, Vec x, Vec r)
	{
		void           *ctx;
		NeutronicsSystem *system ;
		PetscErrorCode ierr = 0;

		// P = gain matrix
		ierr = MatShellGetContext(J, &ctx); CHKERRQ(ierr);
		system = (NeutronicsSystem *)ctx ;

		int n, neig;
		ierr = VecGetSize(x, &n); CHKERRQ(ierr);
		neig = n-1 ;

		// update the solution in NeutronicsSystem so that we know where to 
		// evaluate the nonlinear jacobian when required
		// PetscVector<Real>& true_system_soln = system->get_solution() ;
		// ierr = Wrapper::petsc_copy_vector(x, true_system_soln.vec()) ;CHKERRQ(ierr);
		// Real eigvalue = 0 ;
		// ierr = VecGetValues(x, 1, &neig, &eigvalue) ;CHKERRQ(ierr);
		// system->Keff() = eigvalue ;

		ierr = multiply_jacobian_matrix(J, x, r) ;CHKERRQ(ierr);

		return ierr ;
	}

	// y = J*x for the simple eigenvalue problem
#undef __FUNCT__
#define __FUNCT__ "multiply_jacobian_matrix"
	PetscErrorCode multiply_jacobian_matrix(Mat J, Vec x, Vec r)
	{
		void           *ctx;
		NeutronicsSystem *system ;
		PetscErrorCode ierr = 0;

		// P = gain matrix
		ierr = MatShellGetContext(J, &ctx); CHKERRQ(ierr);
		system = (NeutronicsSystem *)ctx ;

		int n;
		ierr = VecGetSize(x, &n); CHKERRQ(ierr);

		EigenProblemType eigen_type = system->eigen_problem_type() ;

		if (eigen_type == NHEP)
		{
			ierr = Wrapper::multiply_simple_jacobian_matrix(system, x, r) ;
		}
		else if (eigen_type == GNHEP)
		{
			ierr = Wrapper::multiply_generalized_jacobian_matrix(system, x, r) ;
		}
		else
		{
			// Currently, only un-symmetric, simple and generalized eigen problems are handled
			// by the wrapper. All other cases will result in error.
			libmesh_error() ;
		}


		VecAssemblyBegin(r);
		VecAssemblyEnd(r);

		PetscFunctionReturn(ierr);
	}

	// y = J*x for the simple eigenvalue problem
#undef __FUNCT__
#define __FUNCT__ "multiply_simple_jacobian_matrix"
	PetscErrorCode multiply_simple_jacobian_matrix(NeutronicsSystem* system, Vec x, Vec y)
	{
		PetscErrorCode ierr;

		// Jv = [ M*x-Lambda*x-Mu*(Phi) ;
		//			-0.5*(x^T*Phi+Phi^T*x) ] ;
		// Where v = [x; Mu]
		// And Soln = [Phi; Lambda]

		ierr = VecSet(y, 0.0) ;CHKERRQ(ierr);

		int n;
		ierr = VecGetSize(x, &n); CHKERRQ(ierr);

		PetscVector<Real>& pT = dynamic_cast<PetscVector<Real>&>(*system->get_solution().clone().release()) ;
		PetscVector<Real>& pR = dynamic_cast<PetscVector<Real>&>(*system->get_solution().clone().release()) ;

		PetscVector<Real>& pS = system->get_solution () ;	// Phi

		ierr = petsc_copy_vector(x, pT.vec(), false) ;	CHKERRQ(ierr);

		Real jeig = 0.0 ;
		ierr = VecDot(pT.vec(), pS.vec(), &jeig) ;CHKERRQ(ierr);

		jac_context = true ;
		ierr = simple_eigen_operator_private(system, pT.vec(), pR.vec()) ;CHKERRQ(ierr);	// y = M*x
		jac_context = false ;

		Real eigv = system->Keff() ;
		Real neig = 1.0 ;
		//if(KARMA::proc_id() == KARMA::num_procs()-1)
		{
			int nloc ;
			PetscScalar *xarray ;

			VecGetLocalSize(x, &nloc) ;
			VecGetArray(x, &xarray);
			neig = xarray[nloc-1] ;	// Mu
			VecRestoreArray(x,&xarray);
		}

		ierr = VecAXPY(pR.vec(), -eigv, pT.vec()) ;CHKERRQ(ierr);	// y = M*x-Lambda*x 

		ierr = VecAXPY(pR.vec(), -neig, pS.vec()) ;CHKERRQ(ierr);	// pT = pT - Mu*Phi

		ierr = petsc_copy_vector(pR.vec(), y, true) ;CHKERRQ(ierr);	// y = M*x-Lambda*x-Mu*(Phi) 

		//if(KARMA::proc_id() == KARMA::num_procs()-1)
		{
			int nloc ;
			PetscScalar *yarray ;
			VecGetLocalSize(y, &nloc) ;

			VecGetArray(y, &yarray);

			yarray[nloc-1] = jeig ;	// Mu

			VecRestoreArray(y,&yarray);
		}

		// print this information only if solving for preconditioner
		// if (prec_context)
		// {
		// std::cout << "JACOBIAN: keff = " << neig << " and neutronics system keff = " << eigv 
		// << " and eigenvalue jacobian = " << -jeig ;

		// // VecView(x, 	PETSC_VIEWER_STDOUT_WORLD ) ;
		// // std::cin.get() ;
		// }

		Parallel::barrier() ;

		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		PetscFunctionReturn(ierr);
	}

	// y = J*x for the generalized eigenvalue problem
#undef __FUNCT__
#define __FUNCT__ "multiply_generalized_jacobian_matrix"
	PetscErrorCode multiply_generalized_jacobian_matrix(NeutronicsSystem* , Vec , Vec )
	{
		return 0 ;
	}

	// If possible, implement a function to evaluate the diagonal of Shell matrix
#undef __FUNCT__
#define __FUNCT__ "_shell_matrix_get_diagonal"
	PetscErrorCode _shell_matrix_get_diagonal(Mat A, Vec diag)
	{
		void           *ctx;
		NeutronicsSystem *system ;
		PetscErrorCode ierr = 0;

		// A = operator matrix
		ierr = MatShellGetContext(A, &ctx); CHKERRQ(ierr);
		system = (NeutronicsSystem *)ctx ;

		libmesh_assert(system != NULL) ;

		if (system->is_eigen_problem_matrix_free())
		{
			// do nothing for now.
		}
		else
		{
			// PetscVector<Real>& pT1 = system->get_temporary_vector () ;

			// PetscMatrix<Real>* Gain = system->get_gain_matrix() ;

			// // pT1 = P * x ; 
			// ierr = MatMult(Gain->mat(), x, pT1.vec()) ;CHKERRQ(ierr);

			// KSP& loss_solver = system->loss_matrix_solver() ;

			// ierr = KSPSolve(loss_solver, pT1.vec(), y) ;CHKERRQ(ierr);

			// // apply the shift on the operator if any
			// ierr = VecAXPY(y, system->eigen_shift(), x) ;CHKERRQ(ierr);

			// VecAssemblyBegin(y);
			// VecAssemblyEnd(y);
		}

		return 0 ;
	}

	// If possible, apply a shift to the main diagonal of the operator Shell matrix 
	// This is important for shift-invert strategies to work correctly.
#undef __FUNCT__
#define __FUNCT__ "shell_matrix_apply_shift"
	PetscErrorCode shell_matrix_apply_shift(Mat A, PetscScalar shift)
	{
		void           *ctx;
		NeutronicsSystem *system ;
		PetscErrorCode ierr = 0;

		// A = operator matrix
		ierr = MatShellGetContext(A, (void**)&ctx); CHKERRQ(ierr);
		system = (NeutronicsSystem *)ctx ;

		libmesh_assert(system != NULL) ;

		Real& system_shift = system->eigen_shift() ;

		// system_shift += shift ;
		system_shift = shift ;

		// std::cout << "Supplied shift: " << shift << " and updated shift: " << system->eigen_shift() ; 

		PetscFunctionReturn (ierr);
	}

	typedef void (*mat_shell_operations)(void) ;

#undef __FUNCT__
#define __FUNCT__ "shell_matrix_duplicate"
	PetscErrorCode shell_matrix_duplicate (Mat A, MatDuplicateOption op, Mat *M)
	{
		PetscErrorCode ierr = 0;
		void           *ctx;

		PetscFunctionBegin;

		PUSH_LOG("shell_matrix_duplicate", "Eigen Context");

		// A = Jacobian matrix
		ierr = MatShellGetContext(A, &ctx); CHKERRQ(ierr);

		int Nx, Ny, nx, ny ;
		ierr = MatGetLocalSize(A, &Nx, &Ny) ;CHKERRQ(ierr);
		ierr = MatGetLocalSize(A, &nx, &ny) ;CHKERRQ(ierr);

		// For now we care only about 3 shell matrix operations;
		// NOTE: In future if there are more operations that are defined in the neutronics wrapper
		// the function pointers for those operations need to be copied as well.
		mat_shell_operations mult_op, duplicate_op, copy_op, shift_op, zero_op ;

		// Create the shell matrix
		ierr = MatCreateShell(PETSC_COMM_WORLD,nx,ny,Nx,Ny,ctx,M);CHKERRQ(ierr);
		ierr = MatSetFromOptions(*M) ;CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_MULT, (&mult_op)) ;
		ierr = MatShellSetOperation(*M, MATOP_MULT, mult_op);CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_DUPLICATE, (&duplicate_op)) ;
		ierr = MatShellSetOperation(*M, MATOP_DUPLICATE, duplicate_op);CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_SHIFT, &shift_op) ;
		ierr = MatShellSetOperation(*M, MATOP_SHIFT,shift_op);CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_COPY, &copy_op) ;
		ierr = MatShellSetOperation(*M, MATOP_COPY, copy_op);CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_ZERO_ENTRIES, &zero_op) ;
		ierr = MatShellSetOperation(*M, MATOP_ZERO_ENTRIES, zero_op);CHKERRQ(ierr);

		POP_LOG("shell_matrix_duplicate", "Eigen Context");

		PetscFunctionReturn(ierr);
	}

#undef __FUNCT__
#define __FUNCT__ "shell_matrix_copy"
	PetscErrorCode shell_matrix_copy (Mat A, Mat M, MatStructure)
	{
		PetscErrorCode ierr = 0;
		void           *ctx;

		PetscFunctionBegin;

		PUSH_LOG("shell_matrix_copy", "Eigen Context");

		// A = Jacobian matrix
		ierr = MatShellGetContext(A, &ctx); CHKERRQ(ierr);

		// For now we care only about 3 shell matrix operations;
		// NOTE: In future if there are more operations that are defined in the neutronics wrapper
		// the function pointers for those operations need to be copied as well.
		mat_shell_operations mult_op, duplicate_op, copy_op, shift_op, zero_op ;

		ierr = MatShellSetContext(M, ctx) ;CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_MULT, &mult_op) ;
		ierr = MatShellSetOperation(M, MATOP_MULT, mult_op);CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_DUPLICATE, &duplicate_op) ;
		ierr = MatShellSetOperation(M, MATOP_DUPLICATE, duplicate_op);CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_SHIFT, &shift_op) ;
		ierr = MatShellSetOperation(M, MATOP_SHIFT, shift_op);CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_COPY, &copy_op) ;
		ierr = MatShellSetOperation(M, MATOP_COPY, copy_op);CHKERRQ(ierr);

		ierr = MatShellGetOperation(A, MATOP_ZERO_ENTRIES, &zero_op) ;
		ierr = MatShellSetOperation(M, MATOP_ZERO_ENTRIES, zero_op);CHKERRQ(ierr);

		POP_LOG("shell_matrix_copy", "Eigen Context");

		PetscFunctionReturn(ierr);
	}

#undef __FUNCT__
#define __FUNCT__ "precondition_loss_matrix_solve"
	PetscErrorCode precondition_loss_matrix_solve(PC pc, Vec x, Vec y)
	{
		void *ctx ;

		PetscFunctionBegin;
		PetscErrorCode ierr = 0 ;

		ierr = PCShellGetContext(pc,(void**)&ctx);CHKERRQ(ierr);

		NeutronicsSystem* system = (NeutronicsSystem *)ctx ;

		ierr = precondition_loss_matrix(system, x, y) ;CHKERRQ(ierr);

		PetscFunctionReturn(ierr);
	}

	// Shell matrix that represents a preconditioner similar to L
	// Essentially this is a preconditioner that solves L through block-gauss-seidel method
#undef __FUNCT__
#define __FUNCT__ "precondition_loss_matrix"
	PetscErrorCode precondition_loss_matrix(NeutronicsSystem *system, Vec x, Vec y)
	{
		PetscErrorCode ierr = 0;

		PetscFunctionBegin;

		libmesh_assert(system != NULL) ;

		PetscVector<Real> pX(x), pY(y) ;

		if (!prec_created)
		{
			system->recompute_group_blocks() = true ;

			PetscBool point_jacobi = PETSC_FALSE, block_jacobi = PETSC_FALSE ;

			ierr = PetscOptionsHasName("lossoptr_", "-block_jacobi",&block_jacobi);
			if (block_jacobi)
				system->only_block_diagonal_loss_matrix() = true ;
			else
			{
				ierr = PetscOptionsHasName("lossoptr_", "-point_jacobi",&point_jacobi);
				if (point_jacobi)
					system->only_diagonal_loss_matrix() = true ;
			}
		}

		//params.set<bool>("preonly") = prec_prec_context ;

		// block-gauss-seidel preconditioner
		ierr = system->block_gauss_seidel_solve(pX, pY, prec_prec_context) ;  CHKERRQ(ierr);

		// for the sake of it, apply it n times. see what happens ?!
		// {
		// Vec dummy_vec ;
		// ierr = VecDuplicate(x, &dummy_vec) ;  CHKERRQ(ierr);
		// PetscVector<Real> dummy(dummy_vec) ;
		// ierr = system->prec_operator[0]->solve(pY, dummy, params) ;  CHKERRQ(ierr);

		// ierr = system->prec_operator[0]->solve(dummy, pY, params) ;  CHKERRQ(ierr);

		// dummy.clear() ;
		// VecDestroy(&dummy_vec) ;
		// }

		if (!prec_created)
		{
			system->recompute_group_blocks() = false ;
			prec_created = true ;
		}

		PetscFunctionReturn(ierr);
	}

	PetscErrorCode preconditioner_prec_solve(PC pc, Vec x, Vec y)
	{
		void *ctx ;
		PetscErrorCode ierr = 0;

		PetscFunctionBegin;

		ierr = PCShellGetContext(pc,(void**)&ctx);CHKERRQ(ierr);

		NeutronicsSystem* system = (NeutronicsSystem *)ctx ;

		Mat prec_matrix, prec_prec_matrix ;

		int numdofs, locdofs;
		ierr = VecGetSize(x, &numdofs); CHKERRQ(ierr);
		ierr = VecGetLocalSize(x, &locdofs); CHKERRQ(ierr);

		ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)system, &prec_matrix);CHKERRQ(ierr);
		ierr = MatSetFromOptions(prec_matrix) ;CHKERRQ(ierr);

		ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)system, &prec_prec_matrix);CHKERRQ(ierr);
		ierr = MatSetFromOptions(prec_prec_matrix) ;CHKERRQ(ierr);

		ierr = MatShellSetOperation(prec_matrix, MATOP_MULT,(void(*)())Wrapper::preconditioner_multiply);CHKERRQ(ierr);

		KSP ksp_prec ;

		ierr = KSPCreate(PETSC_COMM_WORLD, &ksp_prec) ;CHKERRQ(ierr);
		ierr = KSPSetType(ksp_prec, KSPFGMRES) ;
		KSPAppendOptionsPrefix(ksp_prec,"lossprecp_");
		ierr = KSPSetOperators(ksp_prec, prec_matrix, prec_matrix, SAME_PRECONDITIONER) ; CHKERRQ(ierr);

		PC pc_prec ;
		ierr = KSPGetPC(ksp_prec, &pc_prec) ;CHKERRQ(ierr);
		ierr = PCSetType(pc_prec, PCNONE) ; CHKERRQ(ierr);

		ierr = KSPSetTolerances(ksp_prec, prec_rtol/10, prec_atol, PETSC_DEFAULT, 5) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);
		ierr = KSPSetFromOptions(ksp_prec) ;CHKERRQ(ierr);

		// VecView(x, 	PETSC_VIEWER_STDOUT_WORLD ) ;
		// std::cin.get() ;

		// std::cout << "now printing residual "  ; 
		prec_prec_context = true ;
		ierr = KSPSolve(ksp_prec, x, y) ; CHKERRQ(ierr);

		// VecView(y, 	PETSC_VIEWER_STDOUT_WORLD ) ;
		// std::cin.get() ;
		prec_prec_context = false ;
		ierr = KSPDestroy(&ksp_prec) ;	CHKERRQ(ierr);
		ierr = MatDestroy(&prec_matrix) ;	CHKERRQ(ierr);

		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		PetscFunctionReturn(ierr);
	}


#undef __FUNCT__
#define __FUNCT__ "preconditioner_solve"
	PetscErrorCode preconditioner_solve(PC pc, Vec x, Vec y)
	{
		void *ctx ;
		PetscErrorCode ierr = 0;

		PetscFunctionBegin;

		ierr = PCShellGetContext(pc,(void**)&ctx);CHKERRQ(ierr);

		NeutronicsSystem* system = (NeutronicsSystem *)ctx ;

		Mat prec_matrix, prec_prec_matrix ;

		int numdofs, locdofs;
		ierr = VecGetSize(x, &numdofs); CHKERRQ(ierr);
		ierr = VecGetLocalSize(x, &locdofs); CHKERRQ(ierr);

		ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)system, &prec_matrix);CHKERRQ(ierr);
		ierr = MatSetFromOptions(prec_matrix) ;CHKERRQ(ierr);

		ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)system, &prec_prec_matrix);CHKERRQ(ierr);
		ierr = MatSetFromOptions(prec_prec_matrix) ;CHKERRQ(ierr);

		ierr = MatShellSetOperation(prec_matrix, MATOP_MULT,(void(*)())preconditioner_multiply);CHKERRQ(ierr);

		KSP ksp_prec ;

		ierr = KSPCreate(PETSC_COMM_WORLD, &ksp_prec) ;CHKERRQ(ierr);
		ierr = KSPSetType(ksp_prec, KSPFGMRES) ;
		KSPAppendOptionsPrefix(ksp_prec,"lossprec_");
		ierr = KSPSetOperators(ksp_prec, prec_matrix, prec_matrix, SAME_PRECONDITIONER) ; CHKERRQ(ierr);

		PC pc_prec ;
		ierr = KSPGetPC(ksp_prec, &pc_prec) ;CHKERRQ(ierr);
		ierr = PCSetType(pc_prec, PCSHELL) ; CHKERRQ(ierr);
		ierr = PCShellSetContext(pc_prec, (void*)system) ;
		ierr = PCShellSetApply(pc_prec, (PetscErrorCode (*)(PC,Vec,Vec))preconditioner_prec_solve) ;

		ierr = KSPSetTolerances(ksp_prec, prec_rtol, prec_atol, PETSC_DEFAULT, 20) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);
		ierr = KSPSetFromOptions(ksp_prec) ;CHKERRQ(ierr);

		// VecView(x, 	PETSC_VIEWER_STDOUT_WORLD ) ;
		// std::cin.get() ;

		// std::cout << "now printing residual "  ; 

		ierr = KSPSolve(ksp_prec, x, y) ; CHKERRQ(ierr);

		// VecView(y, 	PETSC_VIEWER_STDOUT_WORLD ) ;
		// std::cin.get() ;

		ierr = KSPDestroy(&ksp_prec) ;	CHKERRQ(ierr);
		ierr = MatDestroy(&prec_matrix) ;	CHKERRQ(ierr);

		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		PetscFunctionReturn(ierr);
	}



#undef __FUNCT__
#define __FUNCT__ "shell_matrix_zero_entries"
	PetscErrorCode shell_matrix_zero_entries (Mat a)
	{
		// do nothing; everything should work just fine !
		// may be set a zeroed out context ??

		PetscFunctionReturn(0);
	}


#undef __FUNCT__
#define __FUNCT__ "preconditioner_multiply"
	PetscErrorCode preconditioner_multiply(Mat A, Vec x, Vec y)
	{ 
		NeutronicsSystem *system ;
		PetscErrorCode ierr=0;
		void           *ctx;

		// A = Jacobian matrix
		ierr = MatShellGetContext(A, &ctx); CHKERRQ(ierr);
		system = (NeutronicsSystem *)ctx ;

		int n;
		ierr = VecGetSize(x, &n); CHKERRQ(ierr);

		EigenProblemType eigen_type = system->eigen_problem_type() ;

		// std::cout << "inside prec * vec .. " ;

		prec_context = true ;
		if (eigen_type == NHEP)
		{
			ierr = Wrapper::multiply_simple_jacobian_matrix(system, x, y) ;
		}
		else if (eigen_type == GNHEP)
		{
			ierr = Wrapper::multiply_generalized_jacobian_matrix(system, x, y) ;
		}
		else
		{
			// Currently, only un-symmetric, simple and generalized eigen problems are handled
			// by the wrapper. All other cases will result in error.
			libmesh_error() ;
		}
		prec_context = false ;

		VecAssemblyBegin(y);
		VecAssemblyEnd(y);

		PetscFunctionReturn(ierr);
	}

	void set_orthonormalize_reference(orthonorm_function user_fn_ref)
	{
		fn_ref = user_fn_ref ;
	}

}


#include "enum_eigen_solver_type.h"

#undef __FUNCT__
#define __FUNCT__ "eigen_residual"
extern "C" PetscErrorCode eigen_residual (SNES snes, Vec x, Vec r, void *ctx)
{
	NeutronicsSystem *system ;
	PetscErrorCode ierr = 0 ;

	system = (NeutronicsSystem *)ctx ;

	// int n, neig;
	// ierr = VecGetSize(x, &n); CHKERRQ(ierr);
	// neig = n-1 ;

	// // update the solution in NeutronicsSystem so that we know where to 
	// // evaluate the nonlinear jacobian when required
	// PetscVector<Real>& true_system_soln = system->get_solution() ;
	// Vec soln ;
	// ierr = SNESGetSolution(snes, &soln) ; CHKERRQ(ierr);
	// ierr = Wrapper::petsc_copy_vector(soln, true_system_soln.vec()) ;CHKERRQ(ierr);
	// Real eigvalue = 0 ;
	// ierr = VecGetValues(x, 1, &neig, &eigvalue) ;CHKERRQ(ierr);
	// system->Keff() = eigvalue ;

	EigenProblemType eigen_type = system->eigen_problem_type() ;

	if (eigen_type == NHEP)
	{
		ierr = Wrapper::simple_eigen_residual(system, x, r) ;
	}
	else if (eigen_type == GNHEP)
	{
		ierr = Wrapper::generalized_eigen_residual(system, x, r) ;
	}
	else
	{
		// Currently, only un-symmetric, simple and generalized eigen problems are handled
		// by the wrapper. All other cases will result in error.
		libmesh_error() ;
	}


	VecAssemblyBegin(r);
	VecAssemblyEnd(r);

	PetscFunctionReturn(ierr);
}

bool is_new_jac_created = false ;
#undef __FUNCT__
#define __FUNCT__ "eigen_jacobian"
extern "C" PetscErrorCode eigen_jacobian (SNES snes, Vec x, Mat *jac, Mat *prec, MatStructure *msflag, void *ctx)
{
	NeutronicsSystem *system ;

	system = (NeutronicsSystem *)ctx ;

	PetscErrorCode ierr = 0 ;

	// if (!is_new_jac_created)
	// {
	// Mat realjac ;
	// ierr = MatComputeExplicitOperator(*jac, &realjac) ;CHKERRQ(ierr);

	// std::cout << "printing realjac matrix " ;
	// PetscViewer view ;
	// PetscViewerCreate(PETSC_COMM_WORLD, &view) ;
	// PetscViewerSetType(view, PETSC_VIEWER_MATLAB) ;
	// ierr = MatCompress(realjac) ;CHKERRQ(ierr);
	// MatView(realjac, 	view ) ;
	// std::cin.get() ;
	// PetscViewerDestroy(&view) ;
	// jac = &realjac ;
	// prec = &realjac ;
	// is_new_jac_created = true ;
	// }

	MatAssemblyBegin(*jac, MAT_FINAL_ASSEMBLY );
	MatAssemblyEnd(*jac, MAT_FINAL_ASSEMBLY );
	MatAssemblyBegin(*prec, MAT_FINAL_ASSEMBLY );
	MatAssemblyEnd(*prec, MAT_FINAL_ASSEMBLY );

	KSP _snes_ksp ;
	SNESGetKSP(snes, &_snes_ksp) ;
	KSPGetTolerances(_snes_ksp, &Wrapper::prec_rtol, &Wrapper::prec_atol, PETSC_NULL, PETSC_NULL) ;

	*msflag = DIFFERENT_NONZERO_PATTERN ;

	// int n, neig;
	// ierr = VecGetSize(x, &n); CHKERRQ(ierr);
	// neig = n-1 ;

	// // update the solution in NeutronicsSystem so that we know where to 
	// // evaluate the nonlinear jacobian when required
	// PetscVector<Real>& true_system_soln = system->get_solution() ;
	// Vec soln ;
	// ierr = SNESGetSolution(snes, &soln) ; CHKERRQ(ierr);
	// ierr = Wrapper::petsc_copy_vector(soln, true_system_soln.vec()) ;CHKERRQ(ierr);
	// Real eigvalue = 0 ;
	// ierr = VecGetValues(soln, 1, &neig, &eigvalue) ;CHKERRQ(ierr);
	// system->Keff() = eigvalue ;

	PetscFunctionReturn(ierr);
}


/*
Monitor - Optional user-defined monitoring routine that views the
current iterate with a petsc viewer. Set by SNESMonitorSet().

Input Parameters:
_snes - the SNES context
its - iteration number
fnorm - 2-norm function value (may be estimated)
ctx - optional user-defined context for private data for the 
monitor routine, as set by SNESMonitorSet()
*/
#undef __FUNCT__
#define __FUNCT__ "snes_monitor"
extern "C" PetscErrorCode snes_monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
{
	PetscErrorCode ierr;

	PetscFunctionBegin;

	int lits = 0 ;
	ierr = SNESGetLinearSolveIterations(snes, &lits) ;CHKERRQ(ierr);

	ierr = PetscPrintf(PETSC_COMM_WORLD,"\nSNES {%D}: \tKSP = %d\t||SNES Residual|| := %G", its, lits, fnorm); CHKERRQ(ierr);


	/*
	NeutronicsSystem* system = (NeutronicsSystem*)ctx ;

	// update the solution in NeutronicsSystem so that we know where to 
	// evaluate the nonlinear jacobian when required
	PetscVector<Real>& true_system_soln = system->get_solution() ;
	Vec soln ;
	ierr = SNESGetSolution(snes, &soln) ; CHKERRQ(ierr);

	int n, neig;
	ierr = VecGetSize(soln, &n); CHKERRQ(ierr);
	neig = n-1 ;

	ierr = Wrapper::petsc_copy_vector(soln, true_system_soln.vec()) ;CHKERRQ(ierr);
	true_system_soln.close();

	if (Wrapper::fn_ref != NULL)
	{
		ierr = (*Wrapper::fn_ref) (true_system_soln.vec(), system->get_temporary_vector().vec()) ;
		CHKERRQ(ierr);
	}
	// true_system_soln.normalize() ;

	Real eigvalue = 0 ;
	ierr = VecGetValues(soln, 1, &neig, &eigvalue) ;CHKERRQ(ierr);
	system->Keff() = eigvalue ;
	*/


	PetscFunctionReturn(0);
}

